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The strongly interacting regime for attractive Bose-Einstein condensates (BECs) tightly confined 
in an extended cylindrical trap is studied. For appropriately prepared, non-collapsing BECs, the 
ensuing dynamics arc found to be governed by the onc-dirncnsional focusing Nonlinear Schrodingcr 
equation (NLS) in the semiclassical (small dispersion) regime. In spite of the rriodulational in- 
stability of this regime, some mathematically rigorous results on the strong asyniptotics of the 
semiclassical limiting solutions were obtained recently. Using these results, "implosion-like" and 
"explosion-like" events are predicted whereby an initial hump focuses into a sharp spike which then 
expands into rapid oscillations. Seemingly related behavior has been observed in three-dimensional 
experiments and models, where a BEC with a sufficient number of atoms undergoes collapse. The 
dynamical regimes studied here, however, arc not predicted to undergo collapse. Instead, distinct, 
ordered structures, appearing after the "implosion", yield interesting new observables that may be 
experimentally accessible. 

PACS numbers: 03.75.Kk, 05.45.Yv 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) with attractive in- 
teractions between atoms have been found to exhibit 
collapse [1] resulting in violent, three-dimensional (3D) 
explosive dynamics [2, 3] and the propagation of quasi- 
one-dimensional (ID), stable bright solitary waves [4, 5]. 
The implosion and subsequent explosion of a 3D BEC, a 
"Bosenova", has been explained theoretically as a blow- 
up singularity of the Gross-Pitaevskii (GP) equation 
when the condensate has a sufHcient number of atoms 
[6-8]. A quasi- ID soliton train was formed by exploiting 
this instability [4]. Other pattern forming instabilities in 
BECs include modulational and transverse instability [9] 
which result in the formation of coherent, localized non- 
linear structures such as solitons and vortices [10]. In this 
work, we explore a new and completely different mech- 
anism leading to violent dynamics and the formation of 
quasi-periodic nonlinear matter wave trains, the semi- 
classical (zero dispersion limit) regime where a quasi-lD 
attractive BEC can "implode" and "explode" yet does 
not undergo 3D collapse. This dynamical regime can be 
accessed with a cigar shaped BEC with negative s-wave 
scattering length of sufRcicntly small magnitude. Recent 
experiments suggest that this may be possible [11]. 

Exact, analytical results for the strong asymptotics of 
the small dispersion limit of the ID focusing (attractive) 
Nonlinear Schrodinger equation (NLS) are used to de- 
scribe the onset of a gradient catastrophe or sharp fo- 
cusing of the condensate density (implosion) followed by 
breaking (explosion) [12, 13]. We emphasize, however, 
that these dynamics do not correspond to collapse of the 
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3D BEC wavefunction. Rather, the resulting dynamics 
reveal two counterpropagating radiative waves, the space 
between them filled by rapidly oscillating quasi-periodic 
(2-phasc) nonlinear matter waves. Further dynamics arc 
determined by the discrete spectrum (solitons) of the 
associated Zakharov-Shabat eigenvalue problem for the 
ID focusing NLS. In the absence of the discrete spec- 
trum, the pure radiation case, the amplitudes of the 
(2-phase) nonlinear matter waves are decaying exponen- 
tially in time [14] . In the presence of the discrete spec- 
trum, (the number of points is inversely proportional to 
the semiclassical parameter) more complicated localized 
structures (2n-phase nonlinear waves) can appear within 
the oscillatory region with n > 1 growing in time. The 
presence of an initial, sufficiently large inward BEC su- 
perfluid velocity (phase gradient) can completely remove 
the localized coherent structures in the oscillatory region 
[12]. 

In this work, we apply the aforementioned rigorous re- 
sults to an attractive quasi-lD BEC in a cigar shaped 
trap. The small dispersion regime can be accessed by 
condensing a sufficiently large number of repulsive atoms 
in a properly sized cigar shaped trap. As has been done 
in the past [2-5, 11], a Feshbach resonance can then be 
applied to tune the sign of the nonlinear term, resulting 
in an attractive mean field interaction. While avoiding 
collapse of the entire condensate, the ID implosion and 
explosion events are shown to be experimentally acces- 
sible; for sufficiently tight radial confinement. Various 
observables resulting from the dynamics of the localized, 
prepared condensate are elucidated including the point 
of gradient catastrophe (breaking point), curves in the 
space-time plane, separating different asymptotic regimes 
(breaking curves), see Fig.l, and the asymptotic struc- 
ture of the BEC density. 

The outline of the paper is as follows. First, we in- 
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troduce the appropriate parameter regime for the ap- 
pUcation of the semiclassical NLS equation to a BEC. 
Following this, we review recent rigorous results on the 
focusing NLS equation with small dispersion and discuss 
their application to single hump initial conditions. 



II. ID BEC AS A SEMICLASSICAL LIMIT OF 
THE FOCUSING NLS 



We now consider the one-dimensional reduction of 
eq. (1) for t > 0. This simplification requires that the 
characteristic energies of the radial excitations are much 
greater than the energies associated with the axial and 
nonlinear excitations [21, 22]. This regime leads to the 
following requirement 
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The temporal evolution of a BEC is governed by the 
GP equation for the condensate wave function ^{r,t), 
given by [15, 16] 



2m 



V^ + V{r,t)+g{t)\^ 



(1) 



where: m is the single atom mass, V = V±{y, z) + V\\{x, t) 
is an external trapping potential with radial V± and axial 
T^l terms, g{t) = 4:Trh^as{t)/m is the nonlinear coefficient 
arising due to two-particle interactions and is character- 
ized by the scattering length Og. We will assume that 



Vu{x,t) 



V|| (x) t < 
t>0 



as{t) 



aP >0 t<0 



(2) 

For t < 0, the BEC is repulsive (positive scattering 
length) and confined in all three spatial dimensions. At 
t = 0, the axial confinement is turned off and the scatter- 
ing length is rapidly switched to a negative value, e.g. by 
a Feshbach resonance [15], so that the condensate be- 
comes attractive. 

We assume that the BEC has been prepared (for t < 0) 
in the ground state of a strongly anisotropic trap 



V±{y,z) = ^muj'}{y^ + z''), 



(3) 



where ujr is the harmonic trap frequency with radial lo- 
calization width ao = \fh[mbj^ The axial portion of the 
potential V|| (a;), for f < 0, confines the EEC's axial extent 
to a width A ^ aq. 

Equation (1) conserves the particle number 



/ 
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A sufficient condition to avoid collapse in a harmonic 
potential is to take A'^ sufficiently small iV < A^cr [17]. 
In two spatial dimensions with a harmonic trap, it has 
been shown that iVcr is related to the loss of stability of 
the nonlinear ground state [7, 18]. Using stability the- 
ory for nonlinear ground states [19], one can compute an 
estimate of N^i- This result has been assumed to hold 
for 3D BECs as well, leading to numerical calculations 
of for several trap configurations [8, 18, 20]. For the 
3D harmonic potential V{x,y,z) = V±{y.z) that we arc 
considering at t > 0, a numerical calculation in [20] de- 
termined 



0.67 



(5) 



where the 0(1) factor C can be estimated through the 
initial density of the BEC (see further discussion in the 
Appendix). If (6) is satisfied, then ^ in eq. (1) for t > 
can be approximated by [21, 22] 




(^{y,z)q{^/A,t/k)e- 



(7) 



where (f>{y,z) = exp[— (j/^ + z^)/(2ao)]/(-v/7rao) is the 2D, 
linear ground state for the transverse harmonic poten- 
tial. Here, the axial variable x has been replaced by ^ for 
convenience. The remaining axial and temporal depen- 
dence is embodied in the function q{x,T) which satisfies 
the NLS equation 
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where, with a slight abuse of notation, we re-introduce 
X = ^/A which is now non-dimensional and set r = t/k. 
The parameters in eq. (8) are 
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We are interested in the semiclassical (small dispersion) 
regime where the semiclassical parameter e satisfies 



< e <Cl. 



(10) 



Conservation of particle number in (4) combined with 
eq. (7) gives 



/ 

Jr 



\q{x, T)\^dx = 1. 



(11) 



In summary, we have derived the NLS equation (8) in 

the small dispersion regime with the assumptions of N < 
(5) and the inequalities (6) and (10). Experimentally, 

all of the parameters ai"\ N, A, and oq can be varied 
so that we expect that these three inequalities can be 
satisfied. By direct calculation, we find that all three 
inequalities (6), (10), and N < Nct can be satisfied by 
choosing, for example, 



N 



ao 



2\ai'''>\ 



<1, 
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2A 



< 1. 



(12) 



The factor C here (and in (6)) is chosen so that 
/r \Q{x,T)\'^dx < C. The inequalities in (12) are just the 
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requirement of cigar shaped initial data with tight radial 
confinement. Using a Feshbach resonance, one can, in 
principle, tune the scattering length ai"' to an arbitrary 
value so that cq. (12) implies that any number of atoms 
is possible. In a recent experiment, the scattering length 
for ''Li has been precisely tuned over seven orders of mag- 
nitude [11]. This demonstrates that the integrable, ID 
NLS equation in the small dispersion limit may be a valid 
model for BEC experiments. 

To illustrate the above discussion, consider the fol- 
lowing data for the bright soliton experiments in [4] 
with ^Li. The parameters involved arc m « 10~^^ kg, 
ai"^ « -1.6 • 10-1° m, ao « 1.6 • 10"^ m, and AT « 3 • 10^. 
With these parameters, A^cr ~ 6700 and the system is 
predicted to undergo collapse, as observed in the ex- 
periment. Nevertheless, assuming an initial axial width 
of A pa 200 /um, which is the approximate experimental 
value, we calculate 

0.011, few 14ms, 7V|4"V^ ~ 0-24. (13) 

Equation (13) ensures the validity of the semiclassical 
regime for at least short times, whereas the quasi-lD as- 
sumption of eq. (13) is, perhaps, on the borderline of 
applicability. 



III. SUMMARY OF PREDICTIONS FROM 
RIGOROUS ASYMPTOTIC ANALYSIS 

Here we bring together and summarize some of the 
main results from the semiclassical, rigorous asymp- 
totic analysis undertaken in [12, 14, 23-27]. Details 
can be found in Sections IV and V. Figure 1 from 
[28] depicts implosion and explosion dynamics gener- 
ated by the single-hump initial condition q{x, 0) = 
exp(— Mncoshx) for the NLS equation (8). Very 
similar dynamics were rigorously derived in [12] for the 
one parameter family of initial data 

g(a;,0,£) = secha;e-5ti"=°^''^, (14) 

where ^ > provides a measure of the phase gradient 
(arg(j)j; or inward supcrfluid velocity in the context of 
BEC [15]. The key features of this study, illustrated by 
Figure 1, can be extended to generic decaying analytical 
one-hump initial data, here are some highlights: 

• The physical x, r plane is subdivided into regions 
where the solution is asymptotically (as £ — )• 0) 
described by modulated 2n-phase nonlinear waves 
(the n = or plane wave approximation corre- 
sponds to the smooth region in Fig. 1, n = 1 to 
the next oscillatory region, etc.). 

• Phase transitions between regions of different be- 
havior are separated by breaking curves in the x, t 
plane that do not depend on e. Equations for the 



breaking curves arc given by (26). A detailed de- 
scription of the transitional behavior at the break- 
ing curve can be found in [37]; 

• The tip of the breaking curve [xo.tq) is the point 
of gradient catastrophe for the plane wave approx- 
imation. Behind this tip, the solution suddenly 
bursts into rapid amplitude oscillations or density 
spikes. Each spike within the vicinity of (a;o,To) 
has the height 'i\q{xQ,Tf))\ and the shape of the ra- 
tional breather solution to the NLS, while the loca- 
tions of the spikes correspond to the poles of the 
special tritronquee solution to the first Painleve 
equation, see [38]. For the family (14), the ex- 
act location of the point of gradient catastrophe 
is {xo,tq) = (0,^) with \q{xQ,TQ)\ = s/ii + 2 so 
that the height of the spikes are 3\//i -|- 2; the slope 
of the breaking curve is -^?==5 . 

• The asymptotic solution for q in the plane wave 
approximation region is completely characterized 
by the implicit formulas in eq. (19) when ji = 2 
and eq. (20) when = 0. In the general case, it is 
determined by eq. (23). 

• When the initial data (14) has ^ < 2, there exist 
©(l/e) points of discrete spectrum (solitons) cen- 
tered at a; = 0. Their interaction leads to the sec- 
ondary break starting at about {x,t) ~ (0, 1.5) in 
Fig. 1, followed by the region of 4-phase wave ap- 
proximation. 

• When the initial data (14) has /x > 2, there are no 
solitons. Therefore, no secondary breaking region 

exists and the solution contains only two regions 
of distinct behavior: smooth and oscillatory. The 
smooth region (n = 0) decays exponentially fast 
in time to zero whereas the background of the os- 
cillatory region decays as 0{t~^/'^) (however, the 
amplitude of the spikes decays exponentially to the 
background with r). A large inward, focusing mo- 
mentum prevents the formation of higher breaking 
regions. 



IV. SEMICLASSICAL LIMIT SOLUTIONS TO 
THE FOCUSING NLS 

In this subsection we discuss some mathematical back- 
ground and recent analytical results related to the semi- 
classical regime of the focusing NLS. The focusing Non- 
linear Schrodinger equation (8) is a universal, basic model 
for self-focusing and self-modulation in that it describes 
the evolution of the envelope of modulated waves in 
generic weakly nonlinear, dispersive systems. It is also 
one of the most celebrated nonlinear integrable equa- 
tions that was first integrated by Zakharov and Shabat 
[29], who used the inverse scattering procedure to de- 
scribe general decaying solutions {\v[n.]^^]^^^ q{x,Q) = 0) 
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FIG. 1: Absolute value \q{x, r, e)| of a solution q{x, t, e) to the 
focusing NLS (8) versus x, r coordinates. Here the initial data 

(15) is given by A{x) = e"'' , S'{x) = — tanhx, and e = 0.02. 
Reprinted from Handbook of Dynamical Systems, Vol 2, D. 
Cai, D. W. McLaughlin, and T. T. R. McLaughhn, The non- 
linear Schrodinger equation as both a PDE and a dynamical 
system, pp 599-675, Copyright (2002) with permission from 
Elsevier. 



in terms of radiation and solitons. Central to their dis- 
covery was the so-called Lax pair which effectively hn- 
earizes the NLS equation. The first equation in the Lax 
pair for the NLS is called the Zakharov-Shabat (ZS) sys- 
tem. It is used to define the correspondence between 
the initial data and the scattering data. Since the scat- 
tering data undergoes an explicit and very simple time 
evolution, the inverse scattering transform (1ST), which 
maps the scattering data back to physical space, is used 
to obtain the evolution of given initial data at any time 
T > 0. 

In the semiclassical limit (e — > 0) , the focusing NLS (8) 
exhibits modulationally unstable behavior (sec Fig. 1), as 
was first shown in [30]. This is in drastic contrast to the 
case of the defocusing (repulsive) NLS equation [28, 31] 
in which the semiclassical theory shows regions of mod- 
ulated periodic or quasi-periodic oscillation. These two 
very different types of behavior can be explained through 
the modiilation eqiiations, which arc elliptic in the fo- 
cusing case and hyperbolic in the defocusing case. The 
corresponding initial value problems are, therefore, ill- 
posed and well-posed respectively. As a result, a plane 
wave with amplitude modulated by A{x) and phase mod- 
ulated by S{x), taken as initial data 



q{x,G,e) = A{x)e 



iS{x)/£ 



(15) 



for the focusing NLS (8), is expected to break immedi- 
ately into disordered oscillations of both the amplitude 
and the phase. However, in the case of analytic initial 
data, the NLS evolution displays some orderly structure 
instead of the disorder suggested by the modulational in- 
stability, see [28, 32, 33]. Throughout this work, we will 
use the abbreviation NLS to mean "focusing Nonlinear 
Schrodinger equation" . 

Figure 1 from [28] depicts the time evolution of a typ- 



ical Gaussian-shaped, symmetric, analytic initial data 
(15). It identifies regions where different types of behav- 
ior (modulated 2n-phase nonlinear waves) of the solution 
q{x, T, e) appear. In particular, consecutive regions with 
n = 0, 2 and, presumably, n = 4 are depicted in Fig. 1. 
These regions are separated by curves in the x, r plane 
that are called breaking curves or nonlinear caustics. The 
location of a breaking curve is defined by A{x) and S{x) 
from (15); it does not depend on e. Within the 2n- 
phase wave approximation region, the strong asymptotics 
of q{x, T, e) can be expressed in terms of Riemann Theta- 
functions (see, for example, [12]), that are defined on the 
genus 2n hyperelliptic Riemann surface TZ{x,t). There- 
fore, the 2n-phase wave approximation region is referred 
to as the genus 2n region. Because of the symmetry 
of the ZS system, TZ{x, r) is Schwarz-symmetrical. The 
surface TZ{x,t) and, more precisely, its 2n -|- 1 complex 
branch points (because of the symmetry, we consider only 
branch points in the upper half-plane), do not depend 
on the semiclassical parameter e. They can be viewed 
as slowly varying functions of the space-time variables 
that describe the wave's parameters, i.e., the macroscopic 
structure of the solution in the vicinity of a given point 

.T, T. 

Equations that define the branch points ofTZ(x, r) are 
known as modulation or Whitham equations. In the case 
n = (genus zero case), TZ{x,t) has only two branch 
points: a = a{x,T) and its complex conjugate a. In 
this case, the Riemann Theta-function expression for 
q{x, T, e) is replaced by 



q{x,T,e) = A{x,T)e 



iS(x,r)/e 



0{e), 



where 



a{x,T) = --S^{x,t) +iA{x,T) 



(16) 



(17) 



and A{x,0) — A(x), S{x,0) — S{x). The genus zero 
region is the first region adjacent to the r = axis where 
the solution (16) has the form of a high frequency C(l/e) 
modulated wave with slowly varying amplitude A{x,t) 
and phase S(x,t). 

For the initial data (14), the corresponding scattering 
data for (8) was calculated explicitly [23] (For general 
initial data of the type (15) sec Sec. V). Using this data, 
the modulated amplitude and phase of (16) and (17) can 
be obtained from the system of transcendental equations 
for a{x, t) = a{x, r) -|- ib{x, r) : 

v/(a - T)2 + 62 + + T)2 + 62 = ^ + 2t62, 
a - r + v/(a - T)2 + bA \a + T+^{a + T)^ + b-^ 

^ L (18) 

where T = 



■ — 1. It is interesting to mention here 

that, in general, the phase gradient S'{x) has a signifi- 
cant impact on the asymptotic behavior of the evolving 
solution. For example [23], in the case < /U < 2, the 
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corresponding ZS eigenvalue problem has O points 
in the discrete spectrum (solitons) located on the vertical 
segment [— T, T\. These solitons are localized at a; = 0, do 
not move, and their interaction can produce quite com- 
plicated coherent structures in the process of evolution. 
In Fig. 1, these structures are seen after the second break 
(around a; = and r > 1.5), which, presumably, corre- 
sponds to the genus 4 region. It is expected that the 
solution undergoes more phase transitions (breaks) for 
larger values of time t that are not visible in Fig. 1. On 
the contrary, in the case > 2, the ZS eigenvalue prob- 
lem does not have any discrete spectrum (solitons). The 
evolution of such initial data produces only two asymp- 
totic regimes (of genera zero and two) , similar to the first 
two asymptotic regions in Fig. 1 (the second asymptotic 
region extends to t = oo, see [12] for the proof). The 
breaking curve separating these two regions has linear 
(slanted) asymptotes as r — >■ oo. In the limit r — >■ c», 
the solution in the genus zero region approaches zero ex- 
ponentially fast, whereas in the genus two region (inside 
the wedge) it decays as 0{t~^) [12]. The high frequency 
amplitude oscillations in this region decay exponentially 
in T, and the solution has the profile of a parabola with 
a maximum at x = and zero values along the break- 
ing curve. Qualitatively, these results mean that a large 
focusing momentum (directed towards the origin) , gener- 
ated by the phase gradient <S"(a;) = — ^ tanhx with /it > 2, 
prevents the formation of the higher genera regions. 

In the borderline case /i = 2, equations (18) have a par- 
ticularly simple solution. Introducing the implicit time 
u = u{x, t) at each point a; G M by t = (u — x) [sinh 2u — 
{u — x)]/(4sinh^ ?i), one can obtain an explicit solution 
of eq. (18) (sec [12]) 

2sinh^u , 2sinhu 

a = . , „ 7 r- , = . , „ 7 r (19) 

smh 2u — [u— X) smh 2u — [u — x) 

for A{x,t) = b and S'{x,t) = —2a that are valid 

throughout the genus zero region. Similar expressions 
with the implicit time u = u{x, r) given by r = 
i a/(u — x) [sinh 2u — {u — x)] coth u and 

9 (u — x) teaih^ u ,n 2tanhM . 

smh 2u — [u — xj smh 2u — [u — x) 

are valid in the case = 0. 

Notice that the amplitude A{x, t) of the solution in 
Fig. 1 at first contracts (focuses) towards the point of 
maximum amplitude, x = 0, and then suddenly bursts 
into rapid (order 1/e) and violent oscillations, transition- 
ing to the genus two regime. This is typical behavior [24] 
for an analytic one-hump initial condition provided that 
S'{x) does not decrease too fast. The very first point of 
this transition, which is the tip-point of the first breaking 
curve (see Fig. 1), is called a point of gradient catastro- 
phe [34]. At the point (xq, tq) of the gradient catastrophe, 
the semiclassical solution (16) of (8) loses its smoothness 
[27], i.e., a^{xn,To) = oc. (cither A^{x,To) or S^^{x,to) 
or both become infinite). The recent results of [38] show 



that each of the spikes seen in Fig. 1 immediately after 
the moment of gradient catastrophe represents a rational 
breather (see [35]) solution to the NLS. The height of each 
spike is exactly three times the value of the amplitude at 
the time of gradient catastrophe, i.e., 3^(a;o, To)+O{e^^^) 
and the location of the spikes are determined by the poles 
of the tritronquee solution of the first Painleve equation. 

Due to the symmetry of the initial data (14), the gra- 
dient catastrophe occurs at xq = 0. In the cases ji = 2 
and /U = 0, the time of the gradient catastrophe tq can be 
calculated as tq = lim„_>.o t(«, 0), where the expressions 
for T = t(u,x) are given above formulae (19) and (20) 
respectively. This yields tq = 1/4 for /i = 2 and tq = 1/2 
for /u = 0. The value of the amplitude bo — A{xo,to) 
at the point of gradient catastrophe can be calculated as 
bo = lim„_>.o &(u, 0), where b = b{u,x) are given in (19) 
and in (20). Thus, bo = 2 for // = 2 and 6o = \/2 for 
fi = 0. In the case /U = 0, numerical solution of (8) with 
e = 1/33 is shown on Fig. 2. 

One possible mechanism for the loss of one- 
dimensionality described by our ID analysis is the gener- 
ation of large axial momentum. This may correspond to 
the breaking down of the energetic assumption of weak 
axial kinetic energy relative to the radial harmonic oscil- 
lator energy (see eqs. (6) and (A. 3)). To illustrate this 
point, we consider the evolution of the Fourier transform 

g(p,T,e) = ^ / y4(a;,r)ei['^(^'^^+^f^lda; (21) 
v27r Js. 

of the solution q(x,T,e) to (8) in the genus zero region 
(where q{x,T,e) has the form (16)). If we define the rel- 
ative momentum P through the standard momentum p 
as P = sp, the standard stationary phase method can be 
applied to (21). Time evolution of \q{P/e,T,e)\ that 

corresponds to the initial data (15) with A{x) = sech x 
and S'{x) = — 2tanha; is shown on Fig. 3 (in this case 
A{x,t) and S{x,t) are given by (17) and (19)). In 
fact, direct calculations show that in the leading e-order 
^lg(P/e,0,e)l = ^ and -1. \q{P/eA/4, e)\ = ^\P\/2 

when \P\ < 2 and ^q{P/e,T,e) = when \P\ > 2 and 
< T < 1/4. Explicit leading order expressions can 
be obtained for any < t < 1/4. They show that the 
portion of atoms in the condensate with high axial mo- 
mentum significantly increases (see Fig. 3) as the point 
of gradient catastrophe at r = 0.25 is approached. 

V. CALCULATION OF SEMICLASSICAL NLS 
SOLUTIONS FOR GENERAL ONE-HUMP 
INITIAL DATA 

Equation (8), the integrable NLS, can be solved by the 
inverse scattering technique. However, the semiclassical 
limit solutions require the semiclassical limit of the scat- 
tering transform. Let E be the curve in the upper half 

plane, defined parametrically by the analytic initial data 
(15) as a{x) = -^^"(a;) + iA{x), x G R. Here A{x) and 
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FIG. 2: Numerical simulation of the solution of (8) with the initial data (14), where jj, — 0, e = 1/33. The time scale t used 
here and our time r are related through r — 2t. This simulation confirm the values to = 0.5 and bo = \/2. It also shows that 
each spike near the point of gradient catastrophe has the height of 36o (with the accuracy (1/33)^/'^ « 0.496) and the shape of 
a rational breather. 



S'{x) — fi±, where fi± are some real numbers, have suf- 
ficient decay as x — >■ ±oo respectively. Let z be a point 
on E. Assuming for simplicity that a{x) is invertible, the 
semiclassical scattering data limit /o(z), z e S, is defined 
[27] through the generalized Abel integral transform as 



z — + ^(z — u)(z — u) x'{u)du 

+ {z-fi+)x{z), (22) 

where x{a) is inverse to a{x) and the integral is taken 
along S. The analytic extension of /o(z) from E to M 
(which can have logarithmic branch cuts) has a meaning 
of the leading order term of ^ie lnro(2:, e) as e — )■ 0, where 
ro{z,e), z S M, is the reflection coefficient of (15). Once 
fo{z) is known, the complex wave parameters encoded 
in the branch points of TZ{x, t) are defined through the 
modulation equations. In particular, in the genus zero 
region, the modulation equation for a{x, r) € is given 
by the system of two real equations [12] 



no 



dC = o, 



Cf(C) 



(23) 



where /(z) = f{z;x,T) = /o(z) — xz — tz"^, 
R{z) = [z — a){z — a) and 7 is an oriented, Schwarz- 



symmetrical contour connecting a and a, such that 
7UM = It defines q(x,T,e) through (16)-(17). Here 
/o(z) is Schwarz symmetrically extended into the lower 
half plane; typically 3/o(z:) has a jump along M. 
We define the function h{z) = h{z]x,T) as 



h{z) 



R{z) 



/(C) 



m L(C-z)i?(C) 



dC - fiz). (24) 



Because of the analyticity of /(z), the particular shape of 
7 is not important. However, it is possible to fix 7 by the 
condition 5/i(z) = on 7. According to the Deift-Zhou 
nonlinear steepest descent method [36], the genus zero 
ansatz (16) approximates the actual solution of the NLS 
(8) with the reflection coefficient ro(z,e) = e~^^°^^^ if 
(see [12]) simultaneously: 

3/i(z; x,t) < on both sides of 7"'"; 
3/i(z;a;,r) > on 7+, (25) 

where j'^ is a contour in the upper half plane con- 
necting a and /i_ and 7"'' = 7nC"'". We have the freedom 
to deform the contours 7,7^ so that the inequalities (25) 
are satisfied along it. The first breaking curve consists 
of points {x, t) where at least one of the inequalities (25) 
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Momentum distribution as t increases in the 
GenusO region 




0.i 



'Mr- 

relative momenUini (P) 



FIG. 3: Evolution of the Fourier transform of the initial data 
(15) with A{x) — sech x and S'{x) = — 2tanhx in the limit 
£ — >■ from r = to the time of gradient catastrophe r — 0.25; 
time t shown on the figure and r are related by t = 2t. The 
vertical axis shows the values of -^l?!, the horizontal axis 
shows the relative momentum P, which is defined as P = 
ep, where p is the standard momentum. The data shown 
here is explicitly calculated from (21) through the standard 
stationary point method. 



turns into equality at some zq. Thus, the equation for 
the first breaking curve can be written as a system of 
three real equations for zq E C and (x, t) e 

9ft.(zo; r) = and hz{zo', x,t) = . (26) 

For the initial data (14) with /i = 2 , the expression 

Va^ + b'^R{z) - a{z - a) + b'^ 



h{z) —z In 



T{z~a)R{z) + z) 
ln[i?(z) — [z — a)], 



ln& 



ITT 

y 



(27) 



was found in [12]. The modulation equations, as well 
as expressions for h(z;x,T) for higher genus regions, can 
be written in explicit determinantal form (see [25, 26]). 
Since these expressions are somewhat involved, they will 
not be given in this paper. 

The point of gradient catastrophe xq , tq is defined by 
the system of equations that consists of (23) and the fol- 
lowing two equations (see [24]): 



no 



cno 

R+iO 



dC = 0. 



This is a system of four real equations for a e 
ixo,To)eR^. 



(28) 



and 



Equations (22)-(28) show that the C'(£)-approximate 
solution in the genus zero region, the point of gradient 
catastrophe, and the breaking curve can be effectively 
calculated from decaying initial data of the form (15). 
Further calculations also reveal the asymptotic structure 
of the solution in the genus two region near the first 
breaking curve and around the point of gradient catas- 
trophe, see [37] and [38] respectively. 



VI. SUGGESTIONS AND CONCLUSIONS 

The semiclassical limit of the focusing ID NLS (8) pro- 
vides a new, mathematically rigorous tool to study the 
modulationally unstable evolution of an attractive ID 
BEC. When the conditions N < Ncr (see eq. (5)) and 
inequahties (6), (lO)-or eq. (12)-are satisfied, an attrac- 
tive BEC in an extended cylindrical trap (cigar shaped 
potential without axial caps) is expected to be governed 
by the focusing NLS in the semiclassical regime. A typ- 
ical evolution for decaying, single-hump, analytic initial 
data q{x,0,e), see (15), is depicted in Fig. 1. Using new 
tools from asymptotic inverse scattering theory, in partic- 
ular, the Deift-Zhou nonlinear steepest descent method, 
outlined above, we were able to calculate the solution 
with 0{e) accuracy in several different regimes. We have 
obtained a number of macroscopic characteristics of the 
evolving condensate. These include (see Fig. 1) the 
space-time location of the point of gradient catastrophe 
and the breaking curves, the slowly modulated ampli- 
tude in the genus zero region, and the envelope of the 
fast amplitude oscillations in the genus two region. 

Attractive ID BEC experiments without axial trap are 
expected to produce two counterpropagating radiative 
waves with the space in between filled with two-phase 
modulated waves and, possibly, with spatially localized 
coherent structures that consist of more complicated non- 
linear waves. These higher genus stationary structures 
are expected to be linked with the discrete spectrum of 
the corresponding ZS system. 

The calculation of the observables, mentioned above, 
is based on the semiclassical limit of the scattering data 
,fo{z), which, in its turn, can be obtained from q{x, 0, e) 
through (22), i.e., through the initial amplitude A{x) and 
the phase S{x). However, accurate measurement of the 
initial phase is often a difficult task. We can turn the 
question around and ask whether the phase S{x) can 
be somehow reconstructed from A{x) and some observ- 
ables. Continuation of this line of argument leads to 
the question of designing some NLS data, initial or scat- 
tering, whose evolution will have certain desired proper- 
ties and/or fit within some required parameters. In light 
of the described above example of the initial amplitude 
A{x) = sech X and the phase gradient S'{x) = — /itanhcc, 
it seems to be especially interesting to observe experi- 
mentally how increasing the phase gradient S'{x) for the 
same amplitude A{x) (by, for example, applying an ex- 
ternal focusing potential) can reduce the number of phase 
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transitions of the evolving condensate. 

Generally speaking, formulae (22)- (26) are valid for a 
large class of analytic initial data, including, for example, 
multi-hump initial data. Further theoretical and experi- 
mental investigations of multi-hump initial densities may 
yield interesting results. 

The authors thank V. Kokoouline for stimulating dis- 
cussions and B. Relethford, who participated in the sum- 
mer 2009 REU DMS 0649159 at the UCF imdcr the su- 
pervision of the first author, for calculating the Fourier 
transform q{p, t, e) and creating Figure 3. 



Appendix: Some energy estimates 

In this appendix, we briefly outline the derivation of 
the quasi-lD criterion (6). 

The conserved energy associated with the 3D GP equa- 
tion (1) for i > (when V{r,t) = Vj_{y,z) and g{t) = 
Awh'^a'f^ /m < 0) is 

£['i/] = + +fni 

«3,l^l^'-*l^ + ^-^^'^')l*l^}^^'+ (A.1) 

2 



approximated by a ID GP equation in the longitudinal 
(axial) direction [21, 22]. For the separated ansatz (14) 
and the ID NLS equation (8) to be valid, we therefore 
assume 



£± » £\\, £± » \£ni\, 
or, since the norm of q is unity (see eq. (11)), 



(A.3) 



^lk.(-,r)||i.(^)«l, ^|k(v 



(A.4) 



Further simplification can be made by using the fact that 
the ID energy for q satisfying the NLS equation (8) is 
conserved 



£ 1 
^iD = y ||fe(-,T)||i2(R) - 2lk(-.^)lli*(K)- (A.5) 

For slowly varying, single-hump (Gaussian type) initial 
data in the semiclassical regime (10) satisfying (11), we 
will have f id < 0. Then the inequalities in (A.3) become 



N\a, 



lk(-,r)||l4(M)-2|£:iD|)«l, 



(«)i 



(A.6) 



< 1. 



Assuming that has the approximate separated form 
given in eq. (7), we formally compute 



Nui±h 2 

^-1 = lk(-,T)||i2(K), 



27r 

47rmA2 



lfe(-,T)||i2(R), 



(A.2) 



£nl = 



2wA ' 



L*(R)' 



where ||/(-)I|lp(m) = (/r l/(a;)|^)^/^ is the standard 
norm. If the characteristic energies of the radial har- 
monic oscillator excitations 6£± > £± are much greater 
then the energies associated with the axial £|| and non- 
linear £ni excitations, the 3D GP equation (1) can be 



If we have 



max(|g(x,T)n<C, Te[0,To], 
then the norm of q is bounded by C because 



(A.7) 



lk(->T)llL4(M) < max(|g(a;,T)| )||g(-,T)||j;,2(R) < C, 

(A.8) 

for T G [0,To] by use of eq. (11). Then the inequalities 
in (A.6) leading to a quasi-lD BEC are satisfied when 
eq. (6) holds. Using (11) and our semiclassical calcula- 
tions, we can estimate C ~ 9|g(xo,To)P for r S [0, Tq] 
with To > To but Tq < ti , the time of the second break. 
For example, the initial data in eq. (14) give C < 9{fi.+2). 
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